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Abstract 

Using a variational approach based on a Lagrangian formulation and Gaussian trial functions, 
we derive a simple dynamical system that captures the main features of the time-dependent 
Schrodinger-Newton equations. With little analytical or numerical effort, the model furnishes 
information on the ground state density and energy eigenvalue, the linear frequencies, as well as 
the nonlinear long-time behaviour. Our results are in good agreement with those obtained through 
analytical estimates or numerical simulations of the full Schrodinger-Newton equations. 
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I. INTRODUCTION 



In recent years, there has been a renewal of interest in the set of nonhnear equations 
known as the Schrodinger-Newton (SN) equations. These consist of the ordinary Schrodinger 
equation 

ih— = A^ + mV(r,t)^, (1) 

at 2m 

where the gravitational potential V{r,t), in the Newtonian approximation, is obtained self- 
consistently from Poisson's equation 

Ay = 47rGm|^|^ (2) 

where m is the mass of the system and G is the gravitational constant. The source term 
in Poisson's equation is provided by a matter density p(r,t) = m|\l/p that is proportional 
to the probability density as given by the wavefunction \I'(r, t). The resulting equations are 
therefore nonlinear. 

SN-type equations have been proposed in various areas of physics and astrophysics. For 
instance, it has been suggested that gravitation, unlike other forces, may not be quantized at 
all fl|. In that case, the stress-energy tensor T^,y in Einstein's equations should be replaced 
by its quantum-mechanical average (T^^)- The SN equations can thus be viewed as the 
nonrelativistic (c — ?■ oo) and Newtonian {G — ?■ 0) limit of the modified Einstein's equations 
G^u = {87fG/c^){Tf^u). More formally, Giulini and GroBardt ^ recently showed that the SN 
equations can be derived in a WKB-like expansion in 1/c from the Einstein-Klein-Gordon 
and Einstein-Dirac system. 



In another context, the SN equations have been proposed as a 



'undamental modification 
4] and Diosi {sj postu- 



of the Schrodinger equation due to gravitational effects. Penrose 
lated that gravity might be at the origin of the spontaneous collapse of the wavefunction 
and proposed the (stationary) SN equations as a possible candidate for an approximate 
description of such gravitationally-induced collapse. 

Finally, the SN equations have been used in an astrophysical context to study self- 
gravitating objects such as boson stars [6] or to describe dark matter by means of a scalar 

field y. 

Whatever their present theoretical status and possible applications, the SN equations 
represent a minimal model in which nonrelativistic quantum mechanics is coupled self- 
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consistently to Newtonian gravity. As such, they are worth investigating in some detail, 
both for their static and their dynamical properties. 

Many theoretical results on the SN equations were obtained in the past, using either 
analytical or numerical approaches 8^]. For instance, the energy eigenvalues (all negative) 
have been determined numerically with good precision |9j and some analytical estimates 
exist on the lower bound for the ground state energy flO|. The linear stability properties of 
the ground state were also investigated jol. 

In the time-dependent and fully nonlinear regime, virtually all results are numerical, with 
few exceptions whose validity is restricted to short time scales An unexpected result was 
published a few years ago by Salzman and Carlip 12|. In numerical simulations of spherically 
symmetric systems, these authors observed that, for masses above a certain critical value, the 
wavefunction "collapsed" at the origin, at least within the accuracy of their simulations 20|. 



The most astonishing feature of these results was that the critical mass was far smaller than 
what could be ex pec ted from simple order-of-magnitude calculations. However, more recent 



calculations 
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13[ disagree with the results of Salzman and Carlip and set the critical 
mass at a value that is several order of magnitudes larger and consistent with analytical 
estimates. 

In this work, we revisit the SN equations using a Lagrangian variational method HQ. 
With this approach, one can arrive at a single ordinary differential equation that describes 
the evolution of the width of the mass density. This method reproduces all the main results 
on the ground state and linear dynamics derived previously. In addition, this approach is 
not restricted to linear theory and can be used to investigate nonlinear oscillations or the 
long-time dynamics. Finally, the mathematical simplicity of the governing equation makes 
it easy to intuit at a glance the salient features of the solutions. 

II. DERIVATION OF THE MODEL 
A. Normalization 

Let us first rewrite the SN equations ([I])-© in dimensionless form, using the analog 
of atomic units for the gravitational interaction. Thus, lengths are measured in units of 
the gravitational "Bohr radius" ac = h^/{Gm?), energy is measured in units of Eg = 
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rrv'G'^ /h^ (the gravitational equivalent of the Hartree), and time in units of Iq = Tl/Eg- To 
ensure conservation of the wavefunction norm, one also needs to normalize ^ to a^?^"^ . The 
dimensionless SN equations then read as 

1 

^^ = -2^^ + nr,t)^, (3) 

A\/ = 47r|^|^ (4) 

with the normalization condition / |\E'p(ir = 1. Notice that Eqs. (I3])-(|1]) are now free of all 
parameters. The evolution of the system is then entirely determined by its initial condition. 
For instance, if the initial condition is spherically symmetric and Gaussian (as will be the 
case in the rest of this paper), the only relevant dimensionless parameter is the width of the 
initial Gaussian measured in units of ac- 



B. Lagrangian approach 

In this section, we will follow the derivation described in Refs. (fsj in the context 
of atomic or condensed matter physics, where the relevant interaction is Coulombian rather 
than gravitational. 

The SN equations ([3])-(jl]) can be written in a hydrodynamical form by using the Madelung 
transformation = y^exp(iS'), where y^p is the amplitude and 5'(r,t) is the phase of the 
wavefunction [igI. The hydrodynamical continuity and momentum equations read as: 

^ + V-(pu) = 0, (5) 

and the velocity is defined as the gradient of the phase, u = VS. 

It can be shown that the above hydrodynamic equations ((5])- ([6]), to geth er with Poisson's 
equation can be derived from the following Lagrangian density C [l5|: 

So far, no approximation was made. The purpose is now to derive a set of evolution equations 
for a small number of macroscopic quantities that characterize the matter density profile. 
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With this aim in mind, let us assume that the system is spherically symmetric and that the 
density profile is Gaussian: 

where r = |r| and R(t) is the time-dependent size of the density. For the above density 
profile, the exact solution of Poisson's equation (jlj) is 

where erf(x) is the error function. In addition, the continuity equation ([5]) is exactly solved 
by the following velocity field: u = {R/R)r, which stems from the phase function S = 
{R/2R)r'^. The dot denotes derivation with respect to time. 

We can now compute the Lagrangian by plugging Eq. ([8]) and the above solutions for V 
and S into Eq. ([7]), and integrating over all space, i.e., L = J Cdr, where the multiplica- 
tive factor was introduced for convenience of notation. The result is 

1-5^5 + 1. (10) 

where C = 2/{3\/2tt). The corresponding equations of motion are obtained from the Euler- 
Lagrange equations 

d dL 
didR 

which yield: 

d^R 

Equation (fT2l) is equivalent to the Hamiltonian equation of motion of a pointlike particle 
evolving in the external potential U{R) = 1/(2/2^) — C/R (see Fig. 1). The first term is 
repulsive and represents kinetic energy due to velocity dispersion (uncertainty principle), 
whereas the second term is attractive and represents self-gravity. 

Note that this result was obtained from a rigorous development based on a Lagrangian 
variational principle. In particular, no assumptions of linear response were made in the 
derivation, so that Eq. ( |T2l) can be used to extract information on the nonlinear regime of 
the time-dependent SN equations. 
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FIG. 1: Radial profile of the pseudo-potential U{R) = 1/{2R^) - C/R. 



III. RESULTS 

The pseudo-potential U{R) is plotted in Fig. 1. It goes to infinity for i? — )■ and goes 
to zero as for R — )■ oo. It crosses the horizontal axis at a point Rq such that U{Rq) = 
and has a single minimum at -Ri, where U'{Ri) = 0. The values of these two points are 
easily determined and yield: 



A. Ground state 

The matter density profile in the ground state is given by Eq. dH]), with R = Ri, corre- 
sponding to the minimum of U. This profile is shown in Fig. 2 (dashed line), together with 
the ground state density obtained from a numerical solution of the stationary SN equations 
(solid line). The agreement is very good. 

Stationary solutions of the SN equations must satisfy the virial theorem, which states 
that the potential energy (in absolute value) is twice the kinetic energy. We can verify that 
this is the case using the wavefunction \I' = y/p from Eq. and the potential of Eq. ([9]). 
The kinetic energy is 




(13) 




(14) 
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FIG. 2: Ground state density: numerically-computed profile (solid line) and Gaussian profile from 
Eq. dlD with R = Ri. 

whereas the potential energy yields: 



It is readily checked that, when R = Ri, \P\ = 2K = l/Svr, so that the virial theorem is 
satisfied. 

The energy eigenvalue of the ground state (lowest energy state)has been computed 
numerically many times and an accepted value is Eq = —0.163 jo], [l^]. In our notation, 
Eq = K + 2P = —3K = —1/2tt ^ —0.159, which is rather close to the numerical value (the 
error is less than 3%). 

Finally, we note that, when R = Rq, one obtains P = —K so that the total energy is 
zero. 

B. Dynamics 

One can linearize the equation of motion f|T2|) around the equilibrium, by writing R = 
Ri + 6R, where 6R(t) is a small perturbation. Substituting into Eq. (fT2|) and taking the 
Fourier transform in the time variable (i.e., assuming 6R{t) ~ e*^*) yields the following 
oscillation frequency: 




(15) 




(16) 
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A perturbation analysis of the full time-dependent SN equations was performed by Har- 
rison et al. jo]. The lowest oscillation frequency that these authors find (see Fig. 2 in Ref. 
9|) is close to ^Harr = 0.035, which is in very good agreement with Eq. (ITB]) (the extra 
factor of 2 comes from the fact that Harrison et al. perturb the wavefunction instead of the 
density I^I/p). 

In order to check this result, we solved the spherically symmetric time-dependent SN 
equations, using a second-order Crank-Nicolson method with centred differences for the 
spatial differentiation. The initial condition is the exact ground state computed numerically 
(solid curve in Fig. 2), to which a very small perturbation was added. The root mean square 



of the radius y (r^it)) is then computed using the standard quantum average. Its frequency 
spectrum is shown in Fig. 3 and displays a clear peak around Q = 0.067, which is again 
very close to Eq. f|T6|l . 




FIG. 3: Frequency spectrum of Vv^WI small oscillations around the ground-state equilibrium, 
obtained from the full time-dependent SN equations. 

For an initial condition that is slightly farther from the exact ground state, we expect 
the SN equations to display some nonlinear effects. This is apparent from Fig. 4, where 
the Gaussian profile given by Eq. with -R(O) = -Ri (dashed line in Fig. 2) was used as 
an initial condition. The time history of \fiyHf)) clearly shows some nonlinear oscillations, 
although their frequency is still close to the linear estimate given by Eq. ( JT6l) . 

Finally, we consider the long-time solutions of Eq. (fT2|) . From the shape of the pseudo- 
potential U{R) (Fig. 1), it is clear that three different regimes are possible for an initial 
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FIG. 4: Evolution of the mean square radius \J ir'^it)) for the full time-dependent SN equations 
with a Gaussian initial condition given by Eq. ([8]) with i2(0) = R\. 

condition i?(0) > 0, i?(0) = 0: 

• If -R(O) < -R05 the total energy is positive, i.e., kinetic energy dominates over grav- 
itational energy. In this case, the wave packet expands indefinitely. Although the 
expansion is slowed down initially by the gravitational attraction, the asymptotic evo- 
lution (t — i- 00) is that of a free particle, i.e., R ^ t. 

• If -R(O) > Ro, the total energy is negative, i.e., gravitational energy dominates over 
kinetic energy. The wave packet oscillates at a nonlinear frequency that can in principle 
be computed from the expression of U (R) (it reduces to the linear frequency Q when 
R{0) ^ 

• If -R(O) = Rq, the total energy is exactly zero. The wave packet still expands, but at 
a rate slower than i? ~ t. The first term on the right-hand side of Eq. ( !T2|) becomes 
negligible for long times. Matching the remaining two terms shows that the expansion 
should go like R ~ t^/^. 

The three regimes described above are neatly reproduced in numerical simulations of Eq. 
( !T2|) . shown in Fig. 5. It is also worth to note that in cosmology these regimes correspond 
respectively to an open, closed, and Einstein-de Sitter universe (which expands as t^/^). 
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FIG. 5: Solutions of the equations of motion (112p for three initial conditions: -R(O) < Rq (red solid 
line), R{0) > Rq (blue line), and R{0) = Rq (black solid line). The dashed straight lines represent 
the curves R ^ t (red) and R ~ t^/^ (black). 

Of course, the full evolution of the wavefunction according to the SN equations can be 
much richer than this simple picture. For large masses, the wavepacket can break down into 
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13|. 



two parts, with some mass being ejected to infinity while the rest remains confined 
This behaviour cannot be captured by our Lagrangian approach, which postulates that the 
density remains close to a Gaussian profile for all times. 



IV. DISCUSSION 



The main interest of the method outlined in this paper is that it relies on a rigorous 
development based on a variational principle, while at the same time yielding results that are 
simple and intuitive. The very shape of the pseudo-potential U{R) (Fig. 1) informs us on the 
type of motions that are to be expected. For instance, it is clear that wavepacket dispersion 
occurs if -R < Rq, whereas it is inhibited if -R > Rq. Further, ii R> Ri 3.7Qfi^ / {Gm?) the 
wavepacket should start to contract right from the beginning of the evolution (for clarity, 
we restore dimensional units in this section). 

Now, we want to compare the above estimations with the numerical results of Giulini et 
al. [ll|, who considered a system of initial size R = 0.707yum (a = 0.5/im in their notation). 
They observed a contracting wavepacket for masses greater than 7 x lO^amu (atomic mass 
units), which is rather close to the value m = 5.74 x lO^amu predicted by our formula 
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Rm^ = 3.76 fi^/G. The results of other simulations 13| are also consistent with these 
findings. 

It is clear that the importance of self-gravitational effects in the SN equations depends 
on both the size R and the mass m of the object under consideration. Therefore, it is useful 
to plot a mass-radius diagram on a log-log scale (Fig. 6), where these two quantities appear 
explicitly. Gravitational effects should play a significant role for objects that fall in the 
region above the curve defined by Rm? = const. = 3.76 fi^/G (solid line). 




10 
m [amu] 



FIG. 6: Mass-radius diagram. Gravitationally induced effects should be important in the region 
above the solid line, which corresponds to the curve Rm^ > 3.76 h^/G. The dashed curve corre- 
sponds to a constant density: m[amu]/i?^ = 5 x lO^^m"'^ (typical solid-state density). 



Experiments aimed at detecting the role of gravity on quantum decoherence will probably 
involve studying the interference fringes of solid-state mesoscopic objects, which should be 
light enough for quantum coherence to be observable but also heavy enough for gravitational 
effects to play a measurable role. Interferometry experiments on small silica spheres 18|] and 
gold clusters [l9| are possible candidates for such studies. 

In order to fix ideas, let us focus on the case of gold or other metal clusters, for which 
the number density is typically rigoid ~ 5 x lO^^m"'^. The dashed line on Fig. 6 represents 
the curve at constant density m[amu]/i?'^ = rigoid- The intersection of the dashed line with 
the the solid line Rm^ = 3.76h^/G yields the minimum mass and radius that gold clusters 
should possess for gravitational effects to play a significant role. This turns out to be of 
the order of a few microns in size and about 5 x 10^ in atomic mass units. The same 
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calculation performed for other metal clusters or the silica spheres mentioned above yields 
similar results. 

The experimental challenge will be to perform quantum interference experiments on such 
massive objects and to control other non-gravitational sources of decoherence. In practice, 
one may perform different experiments for increasing values of the cluster mass, thus moving 
from left to right on the dashed line in Fig. 6. When crossing the solid line, gravitational 
effects should be detected, perhaps as a reduction in the contrast of the interference fringes. 

a. Acknowledgments. F.H. thanks the CNPq (Conselho Nacional de Desenvolvimento 
Cientifico e Tecnologico) for partial financial support. 
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